/*
	Programação Concorrente
	Trabalho 1
	Turma A

	Adilson Rapozo
	Luiz Bertoni
	Rodrigo Costa
*/

#include "../util/util.h"
#include <errno.h>
#include <time.h>
#include <pthread.h>

/*Retorna 1 se ponto esta dentro da circunferencia de raio 1, 0 cc*/
int isInside(long double x, long double y);

/*Thread que executa N iteracoes para algoritmo de Monte Carlo no calculo de PI*/
/*Retorno é a quantidade total de acertos*/
void *monte_carlo(void *arg);

/*Estrutura de parametros da thread*/
/*Contem id da thread e # de iterações que a thread deve executar*/
typedef struct {
	unsigned long id,n_iter;
} parm;

int main(int argc, char* argv[]){

    unsigned long n_total,n_hits,*res,n_process,i;
	mpf_t my_pi;
    clock_t start,end;
	pthread_t *process;
	parm *p;
	FILE *f_print,*f_realpi;

    if (argc != 3){
        printf("usage: %s <# iterações> <# processos>\n",argv[0]);
        exit(1);
    }

    errno = 0;
    if ( argv[1][0]=='-' || (n_total = strtoull(argv[1],NULL,10)) == 0 || errno == ERANGE || errno == EINVAL ){
        printf("Número de iterações não suportado.\n");
        exit(1);
    }

    errno = 0;
	if ( argv[2][0]=='-' || (n_process = strtoull(argv[2],NULL,10)) == 0 || errno == ERANGE || errno == EINVAL ){
		printf("Número de processos não suportado.\n");
		exit(1);
	}

    /*Gera seed para random*/
    srand(time(NULL));


	/*Aloca threads e parametros das threads*/
	process = (pthread_t *)malloc( n_process * sizeof(pthread_t) );
	p = (parm *)malloc( n_process * sizeof(parm) );

    printf("------------------------------\n");

	/*Inicia threads*/
    start = clock();
	for (i=0;i<n_process;i++){
		p[i].id = i;
		p[i].n_iter = (n_total / n_process) + ((i+(n_total%n_process))/n_process); /*espalha uniformemente o # total de iteracoes*/
		pthread_create(&process[i], NULL, monte_carlo, (void*)(p+i));
	}

	/*Aguarda todas threads terminarem e reune os resultados*/
	n_hits = 0;
	for (i=0;i<n_process;i++){
		pthread_join(process[i],(void **)&res); /*res é a qtde de hits retornada pela thread*/
		n_hits = n_hits + (unsigned long)(*res);
		free(res); /*espaço que foi alocado para o retorno da thread*/
	}

	/*Calcula valor obtido de PI*/
    mpf_init2(my_pi,PRECISION_FOR_10M);
    mpf_init_set_ui(my_pi,n_hits);
    mpf_div_ui(my_pi,my_pi,n_total);
    mpf_mul_ui(my_pi,my_pi,4);
    end = clock();

    /*Salva em disco*/
    f_realpi = fopen(REAL_PI,"r");
    if (f_realpi == NULL){
        printf("Arquivo %s não foi aberto.\n",REAL_PI);
        mpf_clear(my_pi);
        free(process);
        free(p);
        exit(1);
    }

    savePI(&f_print,"mc",my_pi,n_total);
    if (f_print == NULL){
        printf("\nArquivo não pode ser salvo.\n");
        mpf_clear(my_pi);
        free(process);
        free(p);
        fclose(f_realpi);
        exit(1);
    }

    /*Saida*/
    printf("------------------------------\n");
    printf("PI = ");
    mpf_out_str(stdout,10,0,my_pi);
    printf("\nWith %lu exact digits calculated right.\n",compare(f_print,f_realpi));
    printf("------------------------------\n");
    printf("Executed in %.2lf seconds.\n",((double)(end-start))/CLOCKS_PER_SEC);

	/*Desaloca variáveis*/
	mpf_clear(my_pi);
    free(process);
    free(p);
    fclose(f_realpi);
    fclose(f_print);
	return 0;
}

/*Thread que executa N iteracoes para algoritmo de Monte Carlo no calculo de PI*/
/*Retorno é a quantidade total de acertos*/
void *monte_carlo(void *arg){
	unsigned long *n_hits,i;
	long double x,y;
	parm *p = (parm *)arg;

    /*numero total de acertos será retornado*/
	n_hits = (unsigned long*)malloc(sizeof(unsigned long));

    /*Inicio das iterações*/
	for(i=0,(*n_hits)=0; i < p->n_iter; i++){
        x = ((long double)rand())/(long double)RAND_MAX;
        y = ((long double)rand())/(long double)RAND_MAX;
        if (isInside(x, y))
            (*n_hits)++;
    }
	printf("Process #%lu finished doing %lu iterations. Hitting %lu times.\n",p->id,p->n_iter,*n_hits);
	return (void *)n_hits;
}

/*Retorna 1 se ponto esta dentro da circunferencia de raio 1, 0 cc*/
int isInside(long double x, long double y){
    if (((x*x) + (y*y)) <= 1)
        return 1;
    return 0;
}

